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Abstract 

The in-plane spatial transport of nonequilibrium excitons in a GaAs quantum 



> 

, well structure has been simulated with the ensemble Monte Carlo method. 

The simulation has been performed for excitons in the presence of residual 
heavy holes including the interparticle Coulomb scatterings, LA phonon scat- 



^ \ terings, and exciton/carrier- interface roughness scatterings. It has been found 

\ that, in contrast to the free electrons/holes system in which the carrier-carrier 

o 



scattering is significant, the interface roughness scattering is the dominant 
process for excitons because of the relatively small scattering rate of exciton- 



H ' carrier and exciton-exciton scatterings. This strongly affects both the spatial 

motion and the energy relaxation of excitons. The spatial and momentum 
distributions of excitons have been simulated up to 500 ps at several exciton 
temperatures and interface roughness parameters. We have found that the 
exciton transport can be regarded as a diffusive motion whose diffusion co- 
efficient varies with time. The diffusivity varies because the average velocity 
of excitons change through the energy transfer between the excitons and the 
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residual heavy hole/lattice system. 
07.05.Tp,02.70.Lq, 73.50.-h, 73.61.Ey, 78.66.-w 



Typeset using REVT^ 



I. INTRODUCTION 



In-plane transport of photogenerated electrons/holes or excitons confined in 2- 
dimensional (2D) space of quantum well (QW) structures has attracted much attention 
as it contains rich physics not accessible in the stationary transport of conduction electrons. 
It is also important from the practical point of view since the spatial resolutions of position 
sensitive optical measurements such as the photoluminescence (PL) via optical microscopes 
and the recently-developed near-field scanning optical microscopes (NSOM)i are determined 
by the diffusion length of photocarriers. 

The measurements of photocarrier transport, which is a transient phenomenon, are dif- 
ficult because it requires high spatial and temporal resolutions. There have been several 
attempts to measure it in QW's of GaAs/AlGaAs with various techniques under different 
experimental conditions. The exciton localization was reported by Hegarty et a/.i using the 
transient grating. Oberhauser et a/.i discussed the scattering mechanism of excitons with 
the same method. The pump-probe technique was utilized by Smith et a/.§ and Yoon et a/.i 
to investigate the excitation density dependence of the spatial motion. Tsen et a/.ifl used 
the time and spatially resolved Raman spectroscopy and discussed the spatial transport as 
well as the energy relaxation. The time-of- flight [tofj method with a photomask was de- 
veloped by Hillmer et a/.i and was applied to study the temperature and excitation energy 
dependence of the transport in QW's under low excitation conditions.iHli The similar tof 
method with an optical fiber was adopted by Akiyama et a/., Ill who reported the transport 
properties in the wide range of excitation density in doped and nondoped QW's. 

In spite of these experimental efforts the physical mechanism of transport is not fully 
understood in the microscopic level so far. This is because the spatial transport of photogen- 
erated carriers is closely related with the longitudinal relaxation of carriers in the momentum 
space (i.e. the energy relaxation). This makes the problem further complicated as compared 
with the transport of conduction electrons for which the momentum space distribution is 
stationary. A distribution function of photocarriers, which are initially generated nonther- 
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mally at some point in the momentum space, changes gradually due to the scatterings of 
carriers, and approaches the thermal equilibrium. At the same time carriers spread spa- 
tially from the point where they are initially generated. Since the scattering processes are 
affected both by the momentum distribution and the spatial distribution (or the density) of 
carriers and the carrier velocities are determined by the momentum distribution, the time 
evolutions of the momentum and the spatial distribution are mutually dependent. Thus it is 
necessary to solve the time evolution of momentum and spatial distribution simultaneously 
to understand the transport properties of photogenerated carriers. Several author^lii re- 
ported the theoretical studies on the spatial transport of excitons in QW's, but they treated 
the momentum distribution in thermal equilibrium, leaving its dynamic evolution out of 
account. 

There are several approaches to the theoretical analysis of the longitudinal momentum 
relaxation (energy relaxation) of photogenerated carriers:lli Analytical estimations based 
on Green's functions, numerical integrations of Boltzmann equations, and numerical simula- 
tions such as ensemble Monte Carlo (EMC) method. The last method is essentially the same 
as solving (nonlinear) Boltzmann equations but, with the advent of fast computers, it has 
become widely used because it is possible to perform quantitative calculations for nonequi- 
librium distributions of carriers and easy to incorporate microscopic scattering mechanisms 
into the model. Another advantage of this method is that the realistic experimental condi- 
tions such as the shape of the sample, the initial, and the boundary conditions, can be used 
in the model calculations, allowing one to compare the theory and the experiment directly. 
(The details of this method is reviewed in Ref. |1^.) 

There are many EMC studies on the nonequilibrium dynamics of photogenerated 
carriers&0 both in bulk and QW structures. But this method simulates the dynamics 
of particles in momentum space and does not yield the information in real space. Another 
class of numerical simulation methods is the Molecular Dynamics (MD) approach, in which 
all particle trajectories are calculated classically in real space. Thus it naturally enables 
one to trace the motion of particles in real space. MD has been applied to the dynamics of 
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nonequilibrium carriers recently by several authors,Ej'Ej however, when the number of par- 
ticipating particles grows, it becomes time-consuming and requires huge computer resources. 
An alternative method is required to practically perform the numerical simulation in the 
long time scale. 

In this paper the spatial transport and the momentum relaxation of excitons is numer- 
ically simulated. The spatial motion is treated by dividing the simulation area by mesh 
(concentric regions in the present case), and the dynamics in the momentum space is calcu- 
lated in each region by EMC method at the density and the momentum distribution averaged 
in the region. The position and the momentum of the particle, and thus the density and the 
momentum distribution in each region are updated at every time step. This approximation 
allows one to trace the spatial distribution of particles at the spatial resolution comparable 
to the region size without the expense of fast computational speed of EMC method. The 
scattering mechanisms as carrier-carrier, exciton-exciton, exciton-carrier, carrier-LA phonon, 
exciton-LA phonon, carrier-interface roughness, and exciton-interface roughness scatterings 
are included in the simulations. 

This paper is organized as follows: In Sec. II the simulation conditions are briefly de- 
scribed, followed by the detailed description the screening and scattering mechanisms, and 
the method of EMC simulation. In Sec. Ill A the comparison of free e/hh and exciton relax- 
ation is discussed and in Sec. Ill B the numerical results of the exciton spatial transport are 
given. They are discussed in comparison with the previous experimental works. Conclusions 
are given in Sec. IV. 

II. MODEL 

A. Simulation Conditions 

The in-plane spatial motion and the energy relaxation of the nonequilibrium excitons are 
calculated in 2D space of nondoped QW structures of GaAs, with the well width Lz = 10 nm. 
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The infinite barrier height is assumed throughout. The simulations have been performed in a 
circular area of radius 10 /im. The initial excitons are spatially distributed in the Gaussian 
shape, C exp ( — 4 In 2 ^^^^^^ ) , with its full width at half maximum (FWHM) 6 yum, and 
their momentum distribution is in the Bose distribution (which is essentially identical to 
the Boltzmann distribution at the density we have used), where the exciton temperature 
Lf,j. is used as a parameter for the initial momentum distributions. The peak areal density 
of excitons is 1.9 x 10^ cm~^ at the center of the simulation area. The residual carriers 
(heavy holes) of uniform areal density 1 x 10^ cm~^ is assumed over the whole simulation 



area as with Ref. |I9|, whose initial momentum distribution is in the Fermi distribution whose 
temperature is the same as the lattice at Tl. [Tl ^ T^x) In Sec. Ill A we also calculate the 
relaxation of free e/hh's to compare with that of excitons. In this case we have used mono- 
energetic initial momentum distributions both for e/hh's and excitons. The simulation is 
performed in the exciton temperatures ^ex) ranging between 10 and 90 K, with the lattice 
temperature {Tl) and the interface roughness parameters varied. The electrons, heavy holes 
and excitons are treated in the parabolic band approximation as it removes the complexity 
of the expressions of carrier-carrier and carrier-LA phonon scattering rates, and contributes 
to the major reduction of computational load. Only the lowest subbands in the F point 
of the conduction and the heavy hole band is considered in the present study, disregarding 
the contribution from the upper subbands and the light hole bands. Since the L-valley is 
located more than 280 meV above, it does not contribute in the energy region of the present 
simulation. If we consider the realistic QW of Alo.3Gao.7As/GaAs with L-^ = 10 nm, the 
energy separations between the lowest hh state and the second hh state or the first light 
hole state is 22 meV and 13 meV, respectively. (The separations are 33 meV and 48 meV, 
respectively, for the well of the infinite barrier height.) Since the average exciton energy 
at the highest exciton temperature (90 K) is 7.7 meV, only a small fraction of excitons in 
the high energy tail of the Bose distribution are involved in the scatterings with the upper 
subbands. Thus we have expected that the upper subbands does not appreciably affect the 
relaxation and transport properties, and we have not included them in the simulations. The 



in-plane {171^^//) and perpendicular {mhh±) (along the growth direction) effective masses in 
the valence band are expressed using Luttinger parameters, 71 and 

rrikh// = (71 + 72)"^ mo, mhh± = (71 - 272)""^ mo 

where tuq is the free electron mass. The physical parameters for GaAs in the present analysis 
are taken from Molenkamp et a/.;iithe effective mass in the conduction band nie = 0.0665mo, 
7i = 6.790, and 72 = 1.924. 



B. Screening 

The proper treatment of the screening is essential in the many-body theory of electrons 
with electron-electron and electron-phonon interactions. There are various approximated 
expressions with the different levels of sophistication.il'0 Among them the random phase 
approximation (RPA) is most commonly used since it predicts important dynamical features 
such as plasmons, and it has analytic expressions in 3D system both for the degenerateii 
and nondegenerateS carrier distributions. But the analytic expression for 2D system is not 
obtained, and we use the static expression in the present study. The original RPA expression 
of the dielectric function, which is valid in 2D for any distributions of carriers, is expressed 



as (see, e.g., Ref. 13^) 



nuj + 10 + J^},_q — 

where is the momentum distribution function of carriers, Ej^ is the energy of the carrier, 
and the summation runs over all possible states of momentum and spin. Vq is the Fourier 
transform of the Coulomb potential in 2D system: 

where is the area of normalization, and Eq is the dielectric constant including the con- 
tributions from the interband electronic transitions and from the phonons. In the static 
[uj —>■ 0) and the long wavelength limit {q —>■ 0) the RPA expression becomes. 



s(q^O, 0) = 1-V,E%^- (2-3) 

k,(T 1 

For 2D, converting the summation over k into an integral over the energy E by utihzing the 
parabohc band dispersion, assuming an isotropic momentum distribution, and integrating 
by parts, Eq. ( p.3| ) reduces to 

e(q - 0, 0) = 1 + E y^i^^f iE = 0), (2.4) 



which is vahd for any isotropic momentum distributions of carriers. Using Eq. (p.2| ) and 
summing up spin indices, one finally obtains 

e{q - 0, 0) = 1 + ^/ = 0) = 1 + -, (2.5) 

Eoh q q 

where the screening wavenumber k is defined by 

«:^^/(E = 0). (2.6) 

This expression for the dielectric function was first adopted by Goodnick et a/.!! and has 
been used for the EMC studies. ii We have used this expression throughout the present study. 
This expression means that only the carriers in the bottom of the band contribute to the 
screening. The screened Coulomb interaction Vg-^'^ is thus described as 

V^ff = Y, = ^ M ^ (2 7) 

^ £(g->0,0) SoL^ \q + ^)' ^ ■ ' 

When the carriers from more than one band have to be considered, they contribute to the 
screening wavenumber independently in the long wavelength limit (g O).0 (Cross terms 
appear when q is finite.) Thus in the present case when the conduction and heavy hole band 
take part in, the screening wavenumber is expressed as 

^ = [fe {E = 0) + hh {E = 0)] (2.8) 

Boh 

where fe (E) and fhh (E) are the momentum distribution functions of electrons and heavy 
holes, respectively. The screening wavenumber is calculated by counting the number of 

8 



electrons and heavy holes at the bottom of each band in every time step of Monte Carlo 
simulations. The strength of Coulomb interaction is updated every time step using this 
screening wavenumber. 

In order to obtain the full carrier dynamics in the relaxation, the dynamic screening 
model is highly desirable. However, the direct calculation of Eq. ( |2.1| ) is computationally 
very heavy and beyond the reach of the present study, (e.g. see Bair and Krusius in Ref. 
351) Instead, the static screening expression of Eq. (^^ ) is used in the following simulations. 

It should be noticed that the expressions for the screening is valid only in the homoge- 
neous system. However, they can be applied to the inhomogeneous system locally when the 
characteristic length of the screening is much smaller than the characteristic length at which 
the distribution functions and the density change.@ In the present study the inverse of the 
screening wavenumber is 100 nm or less, while the characteristic length at which the exciton 
distribution changes is of the order of 1 fim. Thus we can use the locally defined screening 
in the present case. Furthermore, as described below, the excitons do not contribute to the 
screening but only the residual heavy holes, whose distribution is almost homogeneous, are 
responsible for the screening. Thus the effect of inhomogeneity is further reduced. 

We have not included excitons in the screening in the present calculations. The screening 
in the presence of both charged carriers and excitons are discussed in detail by Haug and 
Schmitt- RinkB in the 3D system. In the simplest, case the dielectric function due to excitons 
are given by 

9 

where n and are the exciton density and Bohr radius, respectively. In the present case 
this yields 

£ = eo(l + 5.6 X 10"^). 
On the other hand the dielectric function due to the free carriers is 

e = Soil + 0(1)). 
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Thus the contribution from excitons is much smaller than that from free carriers, and we 
have included only the free carriers in the screening. 



C. Scattering Mechanism 

1. Carrier- Carrier Scattering 
The scattering between charged carriers is discussed in this section. The e-e, hh-hh, and 



e-hh scatterings in 2D using the Born approximation has been discussed in Ref. |T9| in detail 
and only the results relevant to the present study are given here. The total scattering rate 
of an electron of a momentum ki with heavy holes (momentum ^2) is given by 

where is the reduced mass, q is the momentum transfer from the electron to the heavy 
hole, and 9 is the scattering angle in cm. frame. fhh{k2,<y) is the momentum distribution 
function of heavy holes. The screening wavenumber k is given by Eq. (|2.5| ). q is related to 
the initial and the final relative momentum, kr and fc^, and the scattering angle 6 of e-hh 
system by 

q = kr-K, q = 2kr sin(^/2) (2.10) 

where the relative momentum is defined by 

, _ ruhh/fki - me//k2 



rrihh// + rrie// 

Feehh{q) IS the form factor given by 

Feehh{q)= I dze J dzh\Ceiz,)f\Chizh)fe-'^\^^-^'^\, 

—00 —00 

where Ce{ze) and Ch{zh) are the envelope functions in the growth axis for electrons and heavy 
holes, respectively. The form factor is of the order of unity when only the lowest subbands 
of the conduction and the heavy hole bands are considered. Since carrier densities are low 
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in the nondegenerate region, Pauh exclusion in the final state is not considered. The similar 
expression holds for e-e scattering with antiparallel spins. The exchange term is included 
when dealing with e-e scattering with parallel spins. The total scattering rate in this case 
is given byS 



where Q is 



g + K Q + K, 



2 



(2.11) 



Q = 2krCos{9/2). 

The origin of the exchange effect is in the exclusion principle between the identical particles. 
When the two particles are identical, they cannot come close enough. Thus the short range 
interaction does not work between the identical particles. The scale at which the particle 
indistinguishability sets in can be of the order of de Broglie wavelength. In the present case, 
when the relative momentum between two electrons is 0.6 nm~^, the de Broglie wavelength 
is 10 nm. We expect that the exchange effect becomes important when the interaction range 
is equal to or less than the de Broglie length. The interaction range defined by the inverse 
of the screening wavenumber is much longer in our case, around 100 nm. Thus the exchange 
term does not affect much and the interaction strength between the electrons with parallel 
spin is similar to that between the electrons with antiparallel spin. 



2. Electron (Heavy Hole)-LA phonon Scattering 

The LA phonon scattering is treated in the present study since the energy of carriers and 
excitons is less than LO phonon energy and the lattice temperature is low (< 30 K). The 2D 
electrons interact with the phonons of bulk modes through the deformation potential. The 
scattering rate in 2D can be evaluated in parallel with the exciton-LA phonon interaction 
given by Takagahara in detail.il The total LA-phonon (wavenumber Qph) absorption and 
emission rates from the electron with k are given by 
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M/nT(fc) = ^-^2k / de , ^ ,3r (2.12) 



1 D^rr? // 6^ 



+ 1 



(2.13) 



huQph/kT _ ]^ 

where m, p, and Z^e are the velocity of sound, the mass density, and the deformation potential 
for conduction band, respectively. 6 is the electron scattering angle in lab. frame. The simi- 
lar expression holds for heavy hole-LA phonon scattering. Here we assume that the phonon 
energy huQph is much smaller than the electron kinetic energy. The z-component of the 
phonon momentum is always set to zero for simplicity, resulting in the slight overestimation 
of the scattering rate. In the present analysis the physical parameters are taken from Ref. 
13; De = -6.5 eV, Dhh = 3.1 eV, p = 5.3 g cm'^, u = 4.81 x 10^ cm s'^ 

In the present analysis the acoustic phonon distribution is always assumed to be that of 
thermal equilibrium at the lattice temperature, and the effects of nonequilibrium phonons 
emitted from the carriers are not included since the carrier density and the excitation energy 
is low. There are several investigations, theoretically and experimentally,B&^ on the effects 
of nonequilibrium phonons on the dynamics of carriers and excitons. This should be included 
in the case of high density, high excitation energy. 

3. Electron (Heavy Hole) -Interface Roughness Scattering 

The electron (heavy hole)-interface roughness (IFR) scattering in heterostructures was 
given by Ando et a/.0 and applied to the transport phenomena of 2D electron gas in 
QW's.0i3 Here the brief derivation is given for clarity. In the well of infinite barriers, 
the lowest subband energy as a function of well width is given 

When the well width changes from place to place in the 2D plane by ALz^r^/), the subband 
energy changes by 

where r// is the coordinate in 2D plane and me± is the effective mass perpendicular to 2D 
plane. The fluctuation of the well width can be expanded by Fourier series as 
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(2.15) 



'// 



AE{r//) in Eq. ( |2.14 ) can be regarded as a potential for 2D electrons and the e-IFR 
scattering amplitude is given by its matrix element between 2D plane wave states as 



// 



AE(r//) 



k 



II 



(2.16) 



where fc// — h! n = q^^. The Gaussian correlation function of the well- width fluctuation is 
assumed, thus 



(AL,(r//) ■ AL,(r'//)) = exp \{r// - r'//)' /A 



(2.17) 



where A and A are the amplitude and the correlation length of the fluctuation, respectively. 
Using the Fourier expansion Eq. (|2.15|) , the left side reduces to 



The right side can be expanded in Fourier series, then Eq. (|2.17|) reduces to 



// 



'// 



y// 



-g2/A24 



Thus the expansion coefficient is expressed as 

A 



I// 



A 



Using this in the right side of Eq. ( |2.16D , the square of matrix element is written as 



II 



AE{r//) 



k 



II 



2 vr^/i^A^A^ 



-e y/ 



g2 AV4 



(2.18) 



The total scattering rate of the electron with momentum k is calculated via Fermi's golden 
rule as 



2n 



w, 



e-IFR 



>//) 



(2.19) 



where 6 is the scattering angle in Lab. frame. The experimental determination of IFR 
parameters, A and A, is very difficult. The former is the well width fluctuation and usually 
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one or two monolayers is assumed. The latter corresponds to the typical terrace or island 
size in the two dimensional plane. In the present simulations A is 0.283 nm (one monolayer) 
and A is 10 nm unless otherwise stated. In this model, since the infinite barrier height is 
assumed, the IFR scattering is somewhat overestimated compared with the realistic well of 
a finite barrier height. 

4- Exciton- Electron (Heavy Hole) Scattering 

We derive here the exciton-electron (ex-e) scattering to the lowest order. The exciton 
wavefunction in the ideal 2D system (z-dependence ignored) is written as 

where K and R are the center-of-mass momentum and coordinate, respectively. F is the 
wavefunction of relative motion, and c^^ i^Vh) electron (heavy hole) creation operator 

in the Wannier representation. The wavefunction can be rewritten in the Bloch representa- 
tion using the relation 

_ J_ p-ifc-re t 

where N is the number of unit cells in the area A, A = Nvo {vo is the area of the unit cell). 
Converting the summation over lattice sites into the integral by 

E - - / d'r, 

the exciton wavefunction reduces to 

\K) = f{k, k\ K)5j^^j^^yclbl, |0), (2.20) 

where 

/(fe, k\ K)^^J dVF(r)e^('^«^-^)-^. (2.21) 

The mass ratios ctg and ah are defined by 

_ rUe _ rrih 

THe + mh me + THh 
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We evaluate ex-e (or ex-hh) scattering to the lowest order, assuming that there happens 
neither rearrangement of electrons nor the excitations of the internal state of the exciton, 
and ignoring exchange effect. The Coulomb interaction Hamiltonian is written as (spin 
indices are dropped) 

2^1 k+q k q k /c' 

k,k ,q^o 

where is an operator for electrons or heavy holes. V^^^ is the screened Coulomb potential 
of Eq. ( p.7| ), where the screening in only due to the free e/hh's. The scattering amplitude is 
obtained by calculating the matrix element of H' between the initial and the final electron- 
exciton states, \k^K) and \k\K') . Here k {k') and K {K') represent the initial (final) 
electron and exciton momentum, respectively. Two terms remain and the corresponding 
diagrams are shown in (a) and (b) of Fig. 0. The term described by (a) is 



This term can be evaluated by using Eq. ( |2.21|) and converting the summation over momenta 
into the integral. Evaluating two terms (a) and (b), the matrix element is give by 



(fca -q,K + q\H' {k^, K) = ^ / d'r iFfr^^ 



(2.22) 



Using the internal wavefunction of the lowest energy, s-wave exciton; F{r) = 2avoe "''/ v27r, 



where a is the inverse of the Bohr radius, Eq. (|2.22|) reduces to 
V /V^ Aty'^v'^ r r 1 

ZifL^^ / rdrdOe-'"''- e*"'^^'^ - e'^^^^'^ . (2.23) 
j\ 27r J 

This expression can be integrated by using the generation function of the Bessel function, 
and one finally obtains for the matrix element 



{k,-q,K + q\H'\k,,K) = V, 



1 



1 



1 3/2 



1 3/2 



(2.24) 



The total scattering rate of an exciton with electrons (heavy holes) can be calculated with 
the Fermi's golden rule as 
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\ 2 



2-E 



de- 



(q + f^) 



3/2 



1 + 



cteq 
2a 



3/2 



(2.25) 



where is the reduced mass of an exciton and electron (heavy hole), fe,h{k) is the mo- 
mentum distribution function of electrons (heavy holes), and k is the screening wavenumber 
defined in Eq. (p.6| ). The scattering angle 9 is related to the momentum transfer q and the 
relative momentum as g = 2kj./ sin(0/2). The exciton Bohr radius is estimated to 
be 12.5 nm from the variational calculation.Ei 



5. Exciton- Exciton Scattering 



We can derive the exciton-exciton (ex-ex) scattering in the similar manner as ex-e scat- 
tering described in the previous section. The relevant diagrams are shown in (c)-(f) of Fig. 
|l]. The total scattering rate is given as 



W, 



ex—exi-^a) 



16'- 



K 



ex{Ki3 



27r 



d9- 



(q + f^) 



3/2 



1 + 



3/2 



(2.26) 



where is the momentum of the exciton, is the reduced mass of two excitons, and 
fex{K) is the momentum distribution function of excitons. 

Comparing the exciton scattering rates, Eq. ( p.25|) and Eq. ( p.26| ), with the electron 
scattering rate, Eq. ( p.9|) , they show large differences in their magnitude and g-dependence. 
Figure |^ shows the plots of the scattering rates of the three processes vs. q. The e-hh 
scattering rate is forward peaked, and is larger than other two processes by several orders 
of magnitude. The magnitudes of ex-hh and ex-ex scattering rise with q and reach the 
peak at g = 0.15 and 0.20 nm~^, respectively, and are much smaller than e-hh scattering. 
These features of scattering mechanism lead to the striking difference in the relaxation and 
transport between electrons and excitons as described later. 
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6. Exciton-LA Phonon Scattering 



There is a detailed theoretical study of TakagaharaE3 on the exciton-LA phonon inter- 
action and only the results relevant to the present calculations are shown here. The total 
exciton scattering rates for the phonon absorption and emission processes have the similar 
forms as the electron CcLSG cLS, 



4:71 h^upL^ Jo e^"'3ph/fcT _ I 



+ 1 



(2.28) 



where the notations are the same as electron-LA phonon scattering. 



7. Exciton-Interface Roughness Scattering 

The exciton-IFR scattering is deduced from the electron-IFR scattering. Using the ma- 
trix element of Eq. (|2.18|) , the electron-IFR interaction Hamiltonian can be written as 

'rriei, '^//+y// ■"// mhh± '^ii^^ii 



ex-IFR ^1/2^3 ^^^^ fc/z+Q// fe// ^^^^ ^^//' ^ ' 



where and 6^ are operators for an electron and a heavy hole, respectively. The derivation 
is similar to the case of the e-IFR scattering and the total scattering rate is given by 



2lT 



W,^.:pr{K) = ^ ^ dOe-i/'- K (2.30) 

2L° Vme± nihhL' ' 







D. Simulation Method 

The EMC simulations in /c-space have been widely used in the analysis of photogener- 
ated carriers0H^ and reviewed in Ref. |15|, and the detailed technique is not repeated here. 
In order to trace the distributions of the particles in r-space within the frame work of the 
fc-space EMC method, the simulation area is divided into ten concentric circular regions as 
shown in Fig. ^. In each region the average particle density and the momentum distribution 
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function are calculated in every time step and the /c-space EMC simulation is performed, 
independently from other regions, with the scattering parameters corresponding to the aver- 
age density. (Interparticle scatterings are affected directly through their dependence on the 
density and indirectly through the change in the screening wavenumber. Other scatterings 
are also affected indirectly through the change of the momentum distribution of particles.) 
The positions, as well as the momenta, of particles are updated after each time step. A 
particle can move to the neighboring region, thus the particle density, and consequently the 
scattering parameters, change every time step. This technique allows one to trace the time 
evolution of both the spatial and the momentum distribution of particles with relatively 
low computational load. The disadvantage is that, since the continuous density distribution 
(solid line in Fig. |^) is approximated by the step-like one, the difference of the scattering 
rates between the particles near the inner boundary and the outer boundary of the region is 
neglected. Thus the spatial distribution of particles calculated in the present method is only 
valid in the length scale larger than the region size. The range of the Coulomb interaction, 
determined by the inverse of the screening wavenumber is of the order of 100 nm in the 
present simulations, which is much smaller than the region size. 

In the present study ten or twenty simulation areas as mentioned above are prepared and 
simulated simultaneously. The physical information is obtained by taking their ensemble 
average. The number of particles is 78600 when twenty simulation areas are used. The time 
step of simulations are 5 fs for free e/hh's and 12.5 fs for excitons. We have confirmed that 
the time step is much smaller than the scattering interval and that the simulation results 
do not change by further reducing the time step. 



III. RESULTS AND DISCUSSION 
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A. Electron and Exciton Relaxations 



It is expected that the scattering frequency and thus the energy relaxation is different 
between free e/hh's and excitons reflecting the different magnitude of scattering processes. 
To see this, the frequency of each scattering process is counted from to 15 ps both for 
electrons and excitons. For e/hh's the simulation is done in the mono-energetic initial 
momentum distribution, with the excitation energy 10 meV above the band gap (the excess 
energies is 6.34 meV for electrons and 3.66 meV for heavy holes) in the presence of residual 
hh's at the lattice temperature 5 K. The spatial distribution is in the Gaussian function with 
FWHM 6 fim as with excitons. The similar simulation is performed for excitons assuming 
that they were generated mono-energetically at 10 meV. Notice that this is physically not 
realistic because it is not possible to optically generate excitons except at the bottom of the 
band due to the energy-momentum conservation. It is calculated to show the remarkable 
difference between the electron and exciton relaxation due to the difference in the magnitude 
of interparticle scatterings. 

Plotted in Fig. ^ are the numbers of each scattering events per particle per second. For 
electrons e-hh, e-e and e-IFR scatterings are of the similar mag nitude of 1 x 10^^ g-i^ (^jhe 
frequency the e-LA phonon scatterings is several orders of magnitude smaller than those 
plotted here.) The e-hh scattering frequency increases in the initial few ps reflecting the 
rapid energy relaxation of electrons: When the electrons have larger momenta and there 
are many residual hh's with very small average momentum just after photogeneration, e-hh 
scattering is less frequent because it favors the small momentum transfer. 

The exciton scatterings plotted in the lower part of Fig. ^ show different features. (Notice 
that they are plotted in logarithmic scale.) While the ex-IFR scattering is as frequent as 
e-IFR scattering, the interparticle scatterings as ex-ex and ex-hh scatterings are more than 
two orders of magnitude smaller than those of electrons. This is due to the small scattering 
rates of excitons as shown in Fig. |^. Thus the contributions from the interparticle scatterings 
(and LA phonon scatterings which is not shown in Fig. ^ as the figure becomes too busy) is 
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very small. The dominance of IFR scattering in excitons is in agreement with the analytic 
calculations of Basu et a/S 

The difference of the interparticle scattering rates in free e/hh and exciton system has 
significant effect on the time evolutions of the momentum distribution functions. Figure |^ 
shows the distribution functions of the initial 15 ps. The electrons, initially generated at 
0.105 nm^^ (= 6.337 meV), relax rapidly, and approach the thermal distribution described by 
the Boltzmann distribution function (with the electron temperature Te) at 15 ps. This is due 
to the high rate of e-hh and e-e scatterings which are quite efficient in redistributing kinetic 
energies among electrons and heavy holes. It should be noticed that both the electrons and 
the heavy holes are in equilibrium, with the same temperature Tg. (The thermalization 
with the lattice is achieved through the interaction with LA-phonons. But this process is 
quite slow, more than two orders of magnitude lower than e-e scattering. Thus the electron 
temperature is larger than the lattice temperature at this stage.) On the other hand the 
exciton distribution does not change appreciably in the time scale shown here. The ex-IFR 
scattering, which is the dominant process for excitons, is an elastic scattering, and does not 
change the kinetic energy (or the magnitude of momentum) of excitons. Thus the energy 
relaxation is induced only through very slow ex-hh and LA phonon scatterings. Another 
feature attributed to the difference of the scattering processes is that the electron transport 
would depend on the density while the exciton transport would not. 



B. Exciton Transport for 500 ps 

The spatial transport and the energy relaxation of excitons have been calculated up to 
500 ps. Here the initial excitons are in the Bose distribution with the temperature Tex 
while the initial residual heavy holes are in the Fermi distribution at Tl in equilibrium with 
the lattice. Figure |^ shows the time evolution in the momentum distribution functions of 
excitons up to 500 ps with 25 ps step, at Tex = 30 K and the lattice temperature is 5 K. 
Because of the low scattering rates of ex-ex, ex-hh, and ex-LA phonon processes, the energy 
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relaxation is slow, and the excitons do not reach the thermal equilibrium with the lattice 
even at 500 ps. Thus the analytic calculations of the exciton transport which assume the 
momentum distribution in equilibrium with the lattice can underestimate the diffusivity in 
the earlier stage. 

We have obtained the full width at half maximum (FWHM) of the exciton spatial dis- 
tribution by fitting the Gauss function to the distribution calculated from the simulations. 
Figure |^ shows the time evolution of the spatial spread (FWHM) of excitons at several initial 
exciton temperatures {Tex = 10 K - 90 K) from ps to 500 ps when the initial residual hh's 
(and the lattice) are at 5 K. The fitting uncertainty in FWHM is around 0.03 fim. The 
calculations show that the excitons spread faster when the initial exciton temperature is 
higher. This trend can be easily understood when we remember that the dominant scat- 
tering process for excitons is the ex-IFR scattering. The ex-IFR scattering probability (Eq. 
( p.30| )) does not depend on the exciton density and is weakly dependent on the exciton 
momentum. Thus the the mean free path of excitons due to the ex-IFR scattering is similar 
at any place in the simulation area and does not vary significantly during the simulation 
from to 500 ps. In this case, the transport is essentially determined by the velocities of 
excitons, with which they travel during consecutive IFR scatterings. Since the average ve- 
locity is directly related to the exciton temperature, the transport becomes enhanced with 
the temperatures. This also means that, in the short time scale where the exciton temper- 
ature is virtually constant, the transport can be regarded as diffusive motion characterized 
by the constant diffusion coefficient, D = {v)X/2, over whole simulation area. Here (v) is 
an average velocity and A is a mean free path. If excitons were to go through diffusive mo- 
tion with a constant diffusion coefficient and the initial spatial distribution were Gaussian, 
they would spread with time t retaining Gaussian-shaped distribution with its FWHM{t) 
given by y^l6 In 2 x Dt + FWHM{Oy. We have also plotted this curve with the diffusion 
constant 30 cm^ s~^ in Fig. |^ with a broken line for reference. The FWHM with a constant 
diffusion coefficient is almost linear in the scale of Fig. |^, in contrast to the sublinear trend 
of the simulation results. This can be interpreted that the diffusion coefficient decreases 
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dynamically due to the reduction of the average velocity. This is because the initial exci- 
ton temperatures are higher than the initial temperature of residual hh's or the lattice, the 
exciton temperature (and their average velocity) decreases through the ex-hh and ex-LA 
phonon interactions. The diffusion coefficient of the initial stage (0 - 25 ps) is 180 cm^ s~^ 
for Tex = 90 K, while it reduces to 5 cm^ s~^ at 475 - 500 ps. 

In Fig. P we have plotted the spatial spread of excitons when Tl (the initial residual hh 
temperature and the lattice temperature) is varied. In the upper plot, when Tex (30 K) is 
larger than Tl (5 K) (solid curve), the FWHM is sublinear, reflecting the cooling down of 
excitons by the residual hh's and the lattice. When Tex = = 30 K (broken curve), the 
FWHM increases linearly because there is no net energy flow among excitons, residual hh's 
and the lattice. Actually in this case the FWHM evolution can be reproduced by the simple 
diffusion model with a constant diffusivity D = 19 cm^ s^^. On the other hand, in the lower 
plot, when the lattice temperature (30 K) is higher than the initial exciton temperature (10 
K) (broken curve), the FWHM evolution is superlinear because the excitons are heated by 
the residual hh's and the lattice. 

In the calculations above, the IFR parameters were fixed (A = 0.283 nm, A = 10 nm). 
We have done further calculations to investigate how the interface roughness affects the 
transport phenomena. In the present model the interface roughness is characterized by 
the two parameters; A and A. The former designates the typical lateral size of terrace or 
island in the 2D plane of QW's, while the latter indicates the terrace height. There are 
several experimental investigations to observe the interface structures in QW's@'&^l, but 
it is still a controversial problem. We have used several typical values of the parameters in 
the simulations. Figure |^ shows the FWHM of excitons at Tex = 90 K and = 5 K. In 
the upper plot the calculations with A = 5, 10, and 20 nm are compared. (A is fixed at one 
monolayer, 0.283 nm) The diffusion coefficient obtained by fitting to the simple diffusion 
model at the initial stage (0 - 25 ps) is 100, 170, and 250 cm^ s~^ for A = 5, 10, and 20 
nm, respectively. The lower plot shows the A dependence (one monolayer; 0.283 nm, and 
two-monolayer; 0.566 nm) of the transport when A is fixed at 10 nm. The initial diffusion 
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coefficient is 170 an 60 cm^ s~^, respectively. These results show that the transport strongly 
depends on the IFR parameters (the spatial spread is faster when the lateral terrace size is 
larger and the terrace height is lower), and thus the transport properties would vary from 
sample to sample depending on the growth conditions. 

We shall discuss the limitations of the present model with regard to the realistic situation 
in the experiment. The experimentally observed diffusion coefficients vary widely from few 
cm^ s~^ to several hundreds of cm^ s~^ depending on the experimental conditions as sample 
temperatures, excitation energies, carrier densities, and the sample structures. H'ihEl The 
direct comparison of the present simulations with the experimental observation is difficult 
because, in experiments, free electrons and heavy holes are initially photogenerated, then 
they form bound system, excitons. Thus we should have started simulations from free e/hh's 
and handled the exciton formation process properly. However, though there are several 
experimental studies on the exciton formation timeiH^, to the best of our knowledge, the 
theoretical formation model is not available. Here we can compare the simulation results with 
experiments with the similar carrier density conditions only in qualitative manner, assuming 
that if the initial excitation energy is higher, the resultant exciton temperature is also high. 
In Ref. the experimental diffusivity decreases with the excitation energy. The trend 

agrees with the present simulations. As for the sample temperature dependence 
the experimental diffusivity decreases with sample temperatures down to 20 K, which also 
agrees with the simulation results if we assume that the exciton temperature becomes lower 



with the sample temperatures. However, below 20 K, Ref. p!0| , pJ] report the increase of the 
diffusivity, which is not observed in the other experiments. Our simulation results does not 
reproduce the diffusivity enhancement below 20 K. The origin of the enhancement is not 
understood yet. 

For the case of very lower exciton energy, we have to consider the validity of the present 
IFR scattering model. In the ex-IFR scattering an exciton is treated in a plane wave mode, 
extending over all 2D plane. However it is experimentally observed0@ that excitons become 
localized in the potential minima of the interface roughness in 2D plane with the reduction 
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of the in-plane kinetic energy. The locahzed excitons can still migrate in the 2D plane by 
phonon-assisted hopping or variable range hopping (in the very low temperatures). This 
causes the significant slowing down of the diffusivityil. The present model does not include 
this effect. This becomes important when the exciton temperature gets sufficiently low by 
emitting phonons. 



IV. CONCLUSIONS 

The time dependence of the spatial distributions of nonequilibrium excitons was obtained 
using the ensemble Monte Carlo simulations including the interparticle scatterings (ex-ex, 
ex-carrier carrier-carrier), LA phonon emissions/absorptions, and the IFR scatterings. The 
simulations show that the dominant scattering process for excitons is the ex-IFR scattering, 
in contrast to charged carriers for which the carrier-carrier scattering is most significant. 
The difference of the dominant scattering processes affects the energy relaxations: The e/hh 
system approaches the quasi-thermal equilibrium among carriers very rapidly (in few picosec- 
onds) through the efficient energy redistribution caused by the carrier- carrier collisions. In 
contrast, the excitons relax very slowly because their main scattering process, the ex-IFR 
scattering, is elastic and does not change exciton energies. Excitons can exchange energies 
only through the ex-ex, ex-hh and ex-LA phonon scatterings with very low probabilities. 

The transport of excitons is essentially determined by their average velocity because 
the IFR scattering is weakly dependent on the exciton distributions and the effect of other 
scattering processes is negligibly small. Thus the exciton transport can be regarded a diffu- 
sive motion whose diffusion coefficient varies with time. The diffusion coefficient, which is 
proportional to the average velocity, varies through the energy exchange with residual hh's 
and lattice. The spatial transport is strongly dependent on the IFR parameters; the lateral 
terrace size A and the terrace height A. The transport is quenched when the terrace size 
is smaller or the terrace height is higher. This means that the exciton transport properties 
are sensitive to the interface structures, and thus, to the growth conditions of QW's. 
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FIGURES 

FIG. 1. The diagrams calculated in the exciton-electron ((a), (b)) and in the exciton-exciton 

((c)-(f)) scatterings. 

FIG. 2. The comparison of electron-heavy hole, exciton-heavy hole, and exciton-exciton scat- 
tering rates plotted against momentum transfer q. 

FIG. 3. Schematic diagram of the simulation area. The initial excitons are generated in the 
Gaussian shape with FWHM 6 /xm. They gradually spread out to the surrounding area with the 
density distribution described by the solid line. The simulation area (radius 10 /um) is divided into 
concentric circular regions with 1 /xm step, and the average carrier density at each region is used 
in the evaluation of carrier scatterings. 

FIG. 4. The numbers of scattering events per particle per second for electrons (upper) and 
excitons (lower) with the initial excitation energy 10 meV and at 5 K. Notice that in the lower 
case the data are plotted in the logarithmic scale. 

FIG. 5. The time evolution of the momentum distribution functions for (a) electrons and (b) 
excitons, from ps (near side) - 15 ps (far side) with the time step 1 ps. The initial excitation 
energy is 10 meV and the lattice temperature is 5 K. 

FIG. 6. The momentum distribution functions of excitons from ps (near side) to 500 ps (far 
side) with 25 ps step. The initial momentum distribution is in the Bose function with T^.^ = 30 K 
and the lattice temperature at 5 K. 

FIG. 7. The FWHM of the exciton spatial distributions from to 500 ps with the initial exciton 
temperatures changed as a parameter. The lattice temperature is fixed at 5 K. The FWHM is 
obtained by fitting the Gaussian function to the spatial exciton distributions calculated in the 
simulations. The FWHM of the simple diffusion model with D = 30 cm^ s~^ is also plotted for 
reference. 
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FIG. 8. The FWHM of the exciton spatial distributions at the lattice temperature 30 K with 
Tex = 30 K (upper) and 5 K (lower). 

FIG. 9. The FWHM of the exciton spatial distributions at different interface roughness pa- 
rameters: the correlation length A dependence (upper) and the width fluctuation A dependence 
(lower) . 
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